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We propose an efficient method to compute Lyapunov exponents and Lyapunov eigenvectors of 
iong-range interacting many-particie systems, whose dynamics is described by the Viasov equation. 
We show that an expansion of a distribution function using Hermite modes (in momentum variable) 
and Fourier modes (in configuration variable) converges fast if an appropriate scaling parameter is 
introduced and identified with the inverse of the system temperature. As a consequence, dynamics 
and linear stability properties of many-particle states, both in the close-to and in the far-from 
equilibrium regimes can be predicted using a small number of expansion coefficients. As an example 
of a long-range interacting system we investigate stability properties of stationary states of the 
Hamiltonian mean-field model. 
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I. INTRODUCTION 

Neglecting particle correlations, the collective dynam- 
ics of many-particle systems interacting via long-range 
forces obeys the Vlasov equation P, d, y, 0| : 

<9/(x,p,i) <9/(x,p,i) d<p f (x,t) df(x,p,t) 



dt 



dp 



0. 

(1) 

Here /(x, p, t) is the single-particle density function, 
while x and p are, respectively, the particle configura- 
tion and momentum coordinates. The long-range force 
is generated by the self-field 



<p/(x,i) 



J K{x-x')p f {x',t)dx' 



(2) 



where A'(x — x') is the particle interaction potential, and 
p/(x, t) = f /(x, p, t)dp is the particle density. 

There is growing evidence that the evolution of many- 
particle systems with long-range interactions lends itself 
to interpretation in terms of trajectories that trace a fi- 
nite repertoire of collective patterns with abrupt transi- 
tions between them, thus corroborating the claim, stat- 
ing that only a few effectively active degrees of freedom 
partake in dynamics of such systems, despite the huge 
dimensionality of the many-particle state space. 

Presence of transitions indicates also the relevance of 
nonlinear dynamical effects: instabilities, which have the 
capacity to strongly alter the state of a system, even be- 
fore a gradual thermalization of the particles. Analysis 
of a dynamical system in terms of the properties of its in- 
variant structures (equilibria, manifolds, invariant orbits) 
has been established as a fruitful approach to gaining in- 
sights into the nonlinear evolution of a system [f|. In 
particular, the stable equilibria can be identified with rel- 
evant collective patterns, the hyperbolic manifolds of the 
unstable equilibria with pathways of transitions. With 
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the information about the Lyapunov spectrum and Lya- 
punov eigenvectors of such collective patters, control of 
the systems is feasible. 

Effective representation of the equilibria and of their 
local Lyapunov stability properties are one of the main 
challenges addressed in this Paper. Focusing on one- 
dimensional systems represented by density functions 
with Gaussian-like tails in the momentum variable and 
spatially periodic boundary conditions in the configura- 
tion variable, we demonstrate that a Fourier-Hermite ex- 
pansion of f(x,p,t) can reproduce both the Lyapunov 
spectrum and Lyapunov eigenvectors with a small num- 
ber of expansion coefficients. We also show that intro- 
duction of an appropriate scaling parameter, to be iden- 
tified with the inverse of the system temperature, allows 
to improve the convergence of an expansion. 

As a test model, we investigate Lyapunov stability 
properties of stationary states of the Hamiltonian mean- 
field (HMF) model @, 0, I, I M El El El EI EI M ■ 
Prior studies of the HMF model have revealed non-trivial 
collective oscillations [T^] and thermodynamical proper- 
ties that pertain to a large class of long-range interacting 
systems [1 El El El EI El • 

The plan of the Paper is as follows. In Sec. [TT]we intro- 
duce the Fourier-Hermite expansion of a density function 
and discuss the scaling parameter. The spectral equa- 
tion, the Lyapunov spectrum and Lyapunov eigenvectors 
are introduced in Sec. IIIII Relevant facts about the HMF 
model are summarized in Sec. IIV1 Results, obtained from 
the application of our method to the HMF model, are 
presented in Sec. [V] Finally, in Sec. IVI( we draw conclu- 
sions. 



II. THE FOURIER-HERMITE EXPANSION 

In the following, we focus on one dimensional systems. 
Collective properties of many-particle dynamical systems 
are typically more evident when the system states are 
represented in terms of density functions. Further dis- 
cretization of a density function f(x,p,t) in terms of a 



2 



finite set of coefficients a = {a n } provides an alternative 
representation of the system with, hopefully, few cou- 
plings among different a n s. The latter step can be for- 
malized, stating that the density function is a linear map 
H(x,p) of the coefficients a, given by 



f(x,p,t) = U{x,p)a(t) 



(3) 



A Galerkin projection [ljj of a density function f(x,p,t) 
to a set of coefficients is an example of such a map. For 
the Vlasov equation, discretization in momentum vari- 
able commands careful analysis. Therefore we first con- 
sider a partial projection of f(x,p,t): 



f{x,p,t) = ^a n (x,t)yf]^y n (ifip) : 



n>0 



a n {x,t) = 



(4a) 



f(x,p,t)* n (<Jf3p)dp. (4b) 



Here {^f n (y)}n>o are basis functions, and {^ n (y)} n >o 
are basis weights. They constitute complete, normalized, 
adjoint bases, i.e. J 1 $ n (y)^/ m (y)dy = 8 nm . We have in- 
troduced a free scaling function (3, which will be deter- 
mined later. The parameter n n will be fixed by requiring 
that ddn/da n = 0. 

In the following, we make use of the so-called asym- 
metrically weighted Hermite expansion [H, 0j|, [2(| (as 
opposed to the symmetrically weighted Hermite expan- 
sion, discussed in [H[), defined by 

*„(!/) - {2^)- 1 '\-d/dy) n cxp (-y 2 /2) , (5a) 
*"(y) = (n!)" 1 exp (y 2 /2){-d/dy) n exp (-y 2 /2) . (5b) 

Analysis of the Vlasov equation based on Hermite ex- 
pansions in momentum variable has been introduced 
in [iH [2p| . with the emphasis on the accuracy of nu- 
merical solvers. 

Advantages of Hermite polynomials include facts that 
the basis functions decay as ^ n (y) ~ exp (— y 2 /2) for 
large \y\, i.e. in accordance with typical physical bound- 
ary conditions, and that the recursion relations, relevant 
to the Vlasov equation, are simple: 



y* n = A+V n+1 +A n * n - 1 , 



(6a) 
(6b) 



In general, the density pf(x,t) can be expressed as 
Pf(x,t) = a • a, where a = {a n } n >o, a„ = 
J ^> n (y)dy. In the asymmetrically weighted Hermite ex- 
pansion pf(x,t) = ao(x,t). Looking at l[2"]). one can see 
that in this way the couplings between different momen- 
tum modes a n are minimal. This fact motivates the 
choice of the asymmetrically weighted expansion, whose 
recursion coefficients are = 1, A~ = n, B+ = 1, 
B~ = 0, K n = n + 1. 

Using the expansion ((4]), the Vlasov equation (fTJ) is 
cast as 



If 



+ V~ (a)a + V+ (a)a + S~ (a)a + 8+ (a)a = 0, (7) 



where Y^ 1 (a) are defined as 



K(x — x')a ■ a(x' , t)dx' 



(8) 



and 



( a ) = (T)$ T Bt T2 At Tl 8 nT2 , 



BUI 



(9) 



(a n ) - PB 2 A x (a 4 ) 



Here (•} = J -dx. 

The Vlasov dynamics typically proceeds as an inter- 
play between two mechanisms, the advection and the 
convection, and associated with them there are two scales 
of spatial variations. The convection, represented by the 
nonlinear term in |T]), determines the evolution of the 
width of the distribution function in momentum variable. 
It is related to the (time-dependent) system temperature 
T(t): 



m = fp 2 nx,p, t) dxdp. 



(10) 



On this scale, transitions between different macroscopic 
states can be observed. The advection term drives the fil- 
amentation process independently. It is characterized by 
ever-thinning of inhomogeneities (filaments) of /, down 
to the scale of fluctuations in the underlying discrete 
many-particle system. 

In order to represent large-scale variations of / effi- 
ciently, we focus on the former, the temperature scale. 
The parameter (3 in the expansion Q can be tuned to 
maximize the content of the lowest order mode ao(x, t), 
assumed to be nonzero. Using (|4|), one can show that the 
previous condition is equivalent to 



/ 



a^{x, t) dx = , 



[11) 



which in turn corresponds to (3 = 1/T(t). In this way (3 
becomes a dynamical variable, whose evolution is defined 
by enforcing pip . However, for the Lyapunov exponents 
of a stationary state, derivatives of (3, and the two last 
terms in ([7]), can be neglected. 

In the rest of this paper we will study the Lyapunov 
spectrum of smooth distributions with periodic bound- 
ary conditions in x, therefore a Fourier decomposition of 
a n (x,t) is adequate: 



a n {x,t) = 2^ - 



2tt 



exp(- 



(12) 
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With expansions ([4]) and (p~2|) . the system evolution is 
determined by the complex coefficients a = {a mn }, where 
dmn = a -m n- The Vlasov matrix V can be extended to 
include the expansion in x, by noting that the derivative 
d/dx, acting on a mn , is a diagonal operator: d/dx = 
—um. 

For calculations, we use a finite truncation in Np x JVh 
complex coefficients, and let < m < Np and < n < 



ILL LOCAL LYAPUNOV EXPONENTS 



0.15 




0.05 



The local Lyapunov spectrum {(Ji\ and the corre- 
sponding Lyapunov vectors {ipi} of a stationary state ao 
are defined by all solutions (ct,i/>) of the spectral equation 



atp = A(a )ip . 
Here A is the fundamental matrix, defined by 
d 



A(a) 



i— [V"(a)a- 

OSl 



V + (a)a] 



(13) 



(14) 



On a short time scale, the evolution of a perturbed 
state / e (0) = Ila(O) + eU$, (where e is a small arbitrary 
parameter and £ an arbitrary vector) can be written as 
f e (t) w Ila(O) + ell£(i), to the first order in e. The vector 
£ evolves according to the following equation: 



£(t) = Y^l^i ' CX P {WiWi 



(15) 



Here {ipi} is the adjoint basis to {ipi}, satisfying ipitpj = 
Sij. If for some i, A; = Im(— o~i) > 0, the corresponding 
eigenvector ipi is amplified exponentially in time. It is 
therefore referred to as an expanding eigenvector. The 
expansion rate \ is called the local Lyapunov exponent. 
The small number of local Lyapunov exponents is, as 
a rule, a signature of effective low-dimensionality. The 
leading local Lyapunov exponent is defined by 



A = maxlm(— <Xj) 



(16) 



It controls the rate of exponential divergence in time of 
initially proximate states. The divergence is approxi- 
mated by 



Mx,p,t) ^ f (x lP ) + (iP-^e M niP 



(17) 



where ip and tp are the expanding eigenvector and its 
adjoint, respectively, corresponding to the spectral eigen- 
value a with Im(— ct) = A. If all oi are real, the state is 
said to be spectrally stable. 

The rate of exponential divergence between two dis- 
tribution functions is computed by monitoring the evo- 
lution of their difference. This requires the introduction 
of the concept of the distance between two distribution 
functions. We define the distance d(/ e ,/o) as 



<*(/«(*), /o) 




1/2 



(18) 



FIG. 1: Several homogeneous stationary states fo(p;a, (3), 
corresponding to fixed w p = 1.3262. Energies are U a = 0.6, 
U b = 0.72, U c = 0.577 (see also Fig.H) 



IV. THE HMF MODEL 
A. Equations of Motion 

The Hamiltonian mean-field (HMF) model [1, 0, [1] is 
a continuous time model of a one-dimensional "gas" of 

globally coupled particles on a circle, with coordinates 
— 7r < x < 7T, momenta p 6 1R and the "magnetization" 
m(x) = (cosx, smx). The mean magnetization is defined 
by (m) = TV" 1 ^™*, where m, = m(x 2 ). The HMF 
Hamiltonian is 

N 2 N 

ff({ai>M£i) =E| + ^E[ 1 - m .' m ) ]- ( 19 ) 

1=1 — \ 

It determines the particle evolution by 



d 2 , 



df- 



+ (m) A m.; = . 



(20) 



In the limit N — > oo, fixed energy per particle U = H/N, 



(21) 



correlations between particles yield to collective phenom- 
ena. In this, collisionlcss limit, the Vlasov equation for 
the HMF assumes the form 



|Z +p g- m/(f)Am(l) | = 0, (22) 



where f(x,p,t) is the single-particle distribution func- 
tion. The mean-field magnetization is defined by 
m/(t) = J m(x)f(x,p,t)dxdp. 
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FIG. 2: Temporal evolution of a perturbation of the station- 
ary state (a) of Fig. [JJ represented by the distance (|18[) , and 
as a trajectory in the coefficient space in the inset. The thick 
open circle marks the state at t — 20, while the thick crossed 
circle marks the state at t = 50. A fit of the exponential 
regime predicts A = 0.1947 ± 0.0013. The inset shows a pro- 
jection of the trajectory in (ao,2, Re(ai,i), Im(ai,2)). The ini- 
tial perturbation of the stationary state is obtained by Monte 
Carlo integration of the distribution function fo (p) (|23[) with 
one million particles, which is further used as an initial con- 
dition to evolve (l2"Ul l . 



B. Homogeneous Stationary States 

We consider a two-parameter family of homogeneous 
(txif = 0) distributions fo(p;a, (3), given by 



fo(p; a,P) 



1 



1 



2-KWp 1 + exp {(3p 2 /2 - a) 



(23) 



This distribution is a special case of a class of distribu- 
tions, defined by / ~ [1 + exp ((3H — a)]" 1 , shown to 
be important in the statistical treatment of the "violent 
relaxation" processes [2l|. Their relevance to the HMF 
model has been discussed in (l2l |. 

Assuming that fo is normalized, J fo(p) dp = 2ir, 
the parameters a, (3 are related to w v , U by the self- 
consistency conditions [l6[: (3w p = 2n[F_i/ 2 (a)] 2 , U — 
F 1/2 (a)/2(3F_ 1/2 (a) + 1/2. Three examples of f (p; a, (3) 
with different values of U are shown in Fig. [TJ 

The effective temperature /3, determined by (fTTj) , is 



- _ /3F_ 1/2 (a) 
' F 1/2 (a) 



(24) 



where 



1 



r(i + 1) 



exp (t — a) + 1 



are the standard Fermi-Dirac integrals. 

In Ref. [TB], it has been demonstrated that, in the pa- 
rameter plane (w p , U), the stable and the unstable sta- 
tionary states (|23p are separated by the boundary curve 
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FIG. 3: The difference AA between the largest Lyapunov ex- 
ponent ()16|) an d the exact value ()28|) . as a function of the 
truncation order TVh in Hermite polynomials. Here JVf = 5. 
Two cases are shown: the unstable state (a) (dotted line) and 
the stable state (b) (continuous line). The inset shows the 
complete Lyapunov spectrum of the unstable state (a). 



t c (a) = (w P:C (a),U c (a)). Defining Gj(x,a) for a G 
and x e R + by 



Gj(x,a) 



1 



exp(t-a) t j+1 dt 
T(j + l)J [1 + exp (i - a)} 2 t + x 



(25) 



we express i c as 



:(«) 



[^F_ 1/2 (a)G_ 1/2 (0,a)] 1/2 , 

F 1/2 (a)G_ 1/2 (0,a) 1 
4[F_ 1/2 (a)] 2 + 2 



(26) 



In Ref. [16| it has also been demonstrated that for a fixed 
value of w p , U > U m i n (w p ) = 1/2 + w 2 /2A. This defines 
the limiting boundary curve i m i n (w p ) = (w p ,w 2 /2A + 
1/2). Note that approaching £ m ; n corresponds to taking 
the asymptotic limit of a — -> oo of the distribution func- 
tion (|2"3"|) . In this limit fo tends to the "water-bag" dis- 
tribution, defined by f w b{p) = (27rw p ) _1 6(|p| — w p /2), 
where 8(p) is the Heaviside function. Close to this 
boundary, differentiability properties of fo deteriorate, 
and the convergence of the coefficient expansion is slower. 

In the following we will compare our calculations of 
Lyapunov exponents with results on phase transitions for 
the HMF model, discussed in Ref. |12l. Ha|. 

Figure [2] shows the divergence from an equilibrium of 
a perturbed homogeneous stationary state, correspond- 
ing to the curve (a) in Fig. [TJ The evolution has been 
obtained by integrating ([2TJ)l . The initial exponential di- 
vergence is clearly displayed, together with the subse- 
quent saturation to a quasi stationary state, character- 
ized by low- frequency oscillations. The inset of Fig. [2] 
shows the system dynamics, projected in coefficients 
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FIG. 4: Comparison between the exact Lyapunov map and computation using the truncated spectral equation. The left panel 
shows the Lyapunov map (J2SJ . The middle panel shows the results, obtained from the spectral equation (|29[) with Nf = 5, and 
averaging over Nu between 31 and 61. Their difference is shown in the right panel. Special values of parameters, discussed in 
this paper, are indicated by (a), (b), (c). 



(ao,2,Re(ai j i),Im(ai ; 2)). As it can be seen, the stretch 
of the (almost linear) initial trajectory, slightly bending 
towards the attractor just before the reference point at 
t = 20, shows the qualitative features of the long-term 
dynamics. 



V. RESULTS AND DISCUSSION 

Linearization of the Vlasov Equation |T]) around the 
family of distributions (|23p shows that the Lyapunov 
spectrum comprises a continuum of real eigenvalues, as- 
sociated with spectrally stable modes, as well as a fi- 
nite number of imaginary eigenvalues, to be associated 
with unstable collective modes. In the case of homoge- 
neous distributions, Lyapunov exponents of the collective 
modes can be expressed succinctly, using the plasma dis- 
persion function e(a), as an equation e(a) = The 
plasma dispersion function for the HMF model has been 
derived in [8J|: 



e(tr) = 1 + 



-OOP' 



1 df (p) 
dp 



dp, 



(27) 



where /o is given by (|23p . The boundary between stable 
and unstable stationary states is determined by e(0) = 0. 
The result is the curve £ c , given by ip5]l. 

The local (real and positive) Lyapunov exponent A is 
found as a solution of e(i\) = 0. We define this solution 
as the Lyapunov map: 



X(w p , U) 




(28) 



where G 1 (x, a) is the inverse of Gj (x, a) (given by ([2"5])) 
with respect to the variable x, while a = a(w p , U) and 



f3 = /3(w p , U) are expressed in terms of w p and U using 
the self-consistency conditions. 

For non-homogeneous stationary distributions, the dis- 
persion relation does not have a closed form, such as (f27"|) . 
and solving (|13|) is the only direct way to obtain the Lya- 
punov spectrum. Equation (flU)) specialized for the HMF 
model reads 

crip m ,n = m ^ m J~ 1 + (n + l)V> m ,n+l^ 

+ 22 2 i^fi^m-^n-l + a m-k,n-l^k,o] ■ (29) 
fc=±l 

Figure [3] shows the difference between the Lyapunov 
exponent (fT6|) . and the exact value, given by (|28|) . as a 
function of the truncation order iVu in Hermitc polyno- 
mials. The number of Fourier modes is fixed at N-p = 5. 
Two cases are shown, corresponding to the unstable state 

(a) (exact value A a = 0.189549) and to the stable state 

(b) (Ab = 0), in Fig. Q] As it can be seen, TVh > 7 allows 
to predict real Lyapunov spectrum of the stationary state 
(b), and to predict a positive Lyapunov exponent of the 
unstable stationary state (a) qualitatively. Given that 
the sources of errors include both the finite truncation 
error, and errors in the approximation of the coefficients 
(computed using Monte-Carlo integration), the conver- 
gence to the exact value A a can be considered as satisfac- 
tory. The inset of Fig. [3] shows the complete Lyapunov 
spectrum of the unstable state (a), including a single lo- 
cal Lyapunov exponent. All remaining eigenvalues have 
a vanishing imaginary part and tend to cover the real 
axis densely. 

In Fig. [4] we display the Lyapunov map X(w p , U). The 
exact map (|28j) . the approximation computed via (|29|) . 
and their difference are shown in the left, center and right 
panels, respectively. The center panel was computed by 
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FIG. 5: Exponential amplification of a perturbed stationary state (a). The left panel shows the distribution f(t) at t = 20 
(indicated by a thick dot in Fig. [2]) In order to verify (|17(l . the center panel shows the difference between /(20) and the 
initial distribution. The right panel shows the expanding eigenvector, corresponding to A a , computed using (|29|1 with JVf = 5, 
Nu = 11. The phase difference in the configuration coordinate x between the right and middle panels is due to the translation 
invariance of the HMF model. 



averaging the largest Lyapunov exponent over the order 
of truncation, Ah, between A% = 31 and Ah = 61, and 
with the fixed Np = 5. The agreement is fully satisfac- 
tory for the high-energy stable stationary states. Indeed, 
consistently with what is reported in Fig. [3J we have ver- 
ified that for Ah > 7 the computed by (|2Q[) Lyapunov 
spectrum in this region is always real, and therefore pre- 
dictions of spectral stability are unambiguous. Similarly, 
apart from the region close to the intersection between 
l m i n and t c , one can see that also the unstable states are 
predicted unambiguously. The largest Lyapunov expo- 
nent converges to the exact value as the truncation size 
is increased. As for the error in the intersection region, 
analysis, similar to the one shown in Fig. [31 demonstrates 
that close to l c the prediction of A is more sensitive to 
errors in the coefficients (fto)mn- Moreover, deterioration 
of results computed with a fixed truncation order Ah 
close to £ m in can be anticipated because of differentiabil- 
ity properties of /o- 

For comparison, the exact Lyapunov exponent for the 
state (a) A a = 0.1895. Simulation of one million particles 
(see Fig. [5]) gives A = 0.1947. The convergence to A a 
of the calculation based on Fourier-Hermite expansion is 
displayed in Fig. [3J On can see that our method is in 
good agreement with both the A-particle simulation and 
the exact calculation. 

The full solution of (fT3|) provides detailed information 
about local directionality of transitions associated to each 
eigenvalue. In the presence of a single positive local Lya- 
punov exponent, or if one local Lyapunov exponent is 
larger than the remaining ones, this directionality is ex- 
pressed by (jTTJ) . To test accuracy of our method with re- 
spect to calculation of expanding eigenvectors, in Fig. [5] 
we explore the temporal divergence of a perturbed sta- 
tionary state (a). The distribution function, obtained by 
directly integrating (J2DJ) up to t — 20, to near the break- 



down of the regime of exponential divergence (see Fig. [2]), 
is displayed in the left panel of Fig. [6] Subtraction of the 
initial state from the latter is shown in the central panel. 
A direct comparison with the unstable expanding mode, 
calculated using truncated (|29[) with Ap = 11, Ah = 5 
(see the right panel), demonstrates that also Lyapunov 
eigenvectors are predicted accurately. 

As a final step, we tested the ability of the method 
to predict the distribution function after the exponential 
regime, i.e. outside the limit of validity of (|2"9")l . Results 
are summarized in Fig. [S] Considering again the per- 
turbed stationary state (a), the distribution function at 
t = 50 is shown in the left panel of Fig. [5] (indicated by a 
thick crossed circle in Fig. [5]) The effective temperature 
of the distribution is T 1 /3.5. The difference between a 
representation using fixed low-order truncation (Np = 5, 
Ah = 11) is shown in the center panel with (3 = 3.5, and 
in the right panel with f3 = 10. As it can be seen, the 
case in which (3 is tuned at the inverse of the system tem- 
perature results in smaller errors of the approximation. 
The agreement using a detuned j3 can be improved by 
increasing the order of the Hcrmite expansion. 



VI. CONCLUSIONS 

We have proposed an efficient method to compute Lya- 
punov exponents and Lyapunov eigenvectors of collective 
states in long-range interacting many-particle systems, 
whose dynamics is described by the Vlasov equation. The 
method is based on expanding the distribution function 
in a Fourier-Hermite basis, with a scaled momentum vari- 
able. Having tuned the scaling parameter to maximize 
the content of the lowest order coefficient of the expan- 
sion, we have demonstrated that the Lyapunov exponent 
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FIG. 6: The left panel shows the reference distribution f(t) at t = 50, at which time the effective temperature (|10|l is T ~ 1/3.5. 
The center and the right panels show the difference between the reference distribution and the Fourier-Hermite expansion (fj[]), 
with /3 = 3.5 and j3 — 10, respectively. In both cases Nf = 5, Nh = 11. The point of reference t = 50 is also indicated by a 
thick crossed circle in Fig. [5] 



and Lyapunov eigenvectors converge fast to the values 
predicted by direct many-particle simulations and by the 
exact analytical formula. In the last part of the paper, 
the ability of the method to predict the evolution of the 
distribution function in the far-from-equilibrium regime 
has been demonstrated. As an example of a long-range 
interacting system, the Hamiltonian mean-field model 



has been considered. In this context, stability properties 
of a two-parameter family of homogeneous distributions 
over a wide range of dynamical conditions has been inves- 
tigated. Our conclusion is that linear stability properties 
and collective aspects of dynamics of the HMF model 
can be represented and computed using a small number 
of modes, and thus a small numerical effort. 
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